#############################################################################
# Regression of log yields instead of yields
# Basic FE regression
#
#############################################################################

rm(list=ls())
graphics.off()

library(data.table)
library(ggplot2)
library(haven)
library(fixest)
library("RColorBrewer")
library(xtable)
library(Hmisc) # For weighted descriptives

time0=Sys.time()
save_dta = "Data_new"

########################################
# FE estimation
#######################################

# Open data
load("Data_new/FE_results.RData")
save_res = "Results/"
save_results="Results/"

# Generate log yields
summary(dt$yield)
dt[,ly:=log(yield)]
summary(dt[,.(yield,ly)])

# Estimate
fl = feols(ly ~  precr + ddr:factor(fips) | factor(fips) + factor(year)^factor(state),data=dt)
summary(fl)

# Check number of coefficients: 2253 slopes + 1 precr 
nobs(fl)
dim(fl$coeftable)
summary(fl$coeftable)
dt[,ly_hat := fl$fitted.values]
summary(dt[,.(ly,ly_hat)])

# Marginal effects in logs (mgl) and in levels (mg)
source("Code/0. CODE Auxiliary.R")
dt_fel = slopes_dt(f=fl,names_var="")
dt_fel = merge(dt[,.(fips,year,yield,ly,ly_hat)],dt_fel[,.(fips,coef)],by="fips",all.y=TRUE)
setnames(dt_fel,"coef","beta")
dt_fel[, mg_log:= beta*exp(ly_hat+fl$sigma2/2)]
dim(dt_fel)
summary(dt_fel)

data_all = dt_res[,.(fips,year,mg3,corn_area)]
dim(data_all)
data_all = merge(data_all,dt_fel[,.(fips,year,mg_log)],by=c("fips","year"),all=TRUE)
dim(data_all)
summary(data_all)
setnames(data_all,"mg3","fe_het")

time1=Sys.time()

########################################
# Describe results
########################################

# Density
g = ggplot(data=data_all) + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18), axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = 'right',legend.text = element_text(size=14)) +
  scale_color_discrete(name= "Groups") +
  xlab("Temperature effects") + ylab("Density") 
yl=c(0,0.16)
xl=c(-30,10)
g + geom_density(aes(x=fe_het,weight=corn_area),adjust=1.5,linetype="dashed") +
  geom_density(aes(x=mg_log,weight=corn_area),adjust=1.5,linetype="solid") +
  xlim(xl) + ylim(yl)
ggsave(paste0(save_res,"Density_fe_logs.pdf"),width=6.3,height=3.9,units="in")

# Map
source("Code/0. CODE Auxiliary.R")
load("Data_ori/Data_us_map.RData")

dt_plot = data_all[is.na(fips)==FALSE,lapply(.SD,wtd.mean,weights=corn_area),by="fips",.SDcols=c("fe_het","mg_log")]
dim(dt_plot)
summary(dt_plot)
qq= quantile(dt_plot$mg_log,probs=seq(0.05,0.95,0.05))
qq
bmin= qq[2]
bmax= qq[18]

# FE long panel
plot_map(data=dt_plot[,.(fips,mg_log)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_mg_fe_log",palette="OrRd",rev=1,nc=9)


